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ABSTRACT 

This  paper  describes  a  successive  overrelaxation  (SOR)  method  for  computing  a  class 

,  \ 

of  bivariate  C*  piecewise  cubic  polynomial  interpolants.  Given  a  collection  of  points  in  R* 
together  with  a  triangulation  of  those  points,  the  schemes  described  require  only  the  values 

i 

of  the  function  to  be  interpolated  at  the  given  points.  The  result  is  a  C J  interpolant  which 
is  a  cubic  polynomial  over  each  of  the  triangles  in  the  triangulation.  Numerical  results  are 
presented  for  a  typical  scheme  in  this  class. 
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SIGNIFICANCE  AND  EXPLANATION 


This  paper  describes  a  class  of  bivariate  C1  piecewise  cubic  polynomi-'  interpolation 
schemes  for  scattered  data.  Computation  of  this  interpolant  requires  the  solution  of  a 
linear  system  of  equations  which  is  extremely  sparse.  Since  the  triangulation  of  scattered 
data  may  allow  for  any  number  of  triangles  to  meet  at  a  vertex,  the  matrix  corresponding 
to  this  linear  system  need  not  have  a  nice  structure.  This  means  that  the  ordering  of  the 
vertices  of  the  triangulation  may  be  of  critical  importance  whenever  a  direct  method  of 
solving  the  linear  system  is  to  be  employed.  In  this  paper,  an  iterative  method  of  solving 
this  system  is  proposed  in  order  to  eliminate  this  difficulty.  Numerical  evidence  suggests 
that  this  iterative  scheme  is  an  extremely  efficient  and  effective  way  of  solving  the  linear 
system. 
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The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary  lies 
with  MRC,  and  not  with  the  author  of  this  report. 


AN  ITERATIVE  METHOD  FOR  COMPUTING  BIVARIATE  C1 
PIECEWISE  CUBIC  POLYNOMIAL  INTERPOLANTS 


Thomas  A.  Grandine1,2 


Interpolation  of  bivariate  scattered  data  by  smooth  piecewise  polynomials  is  a  problem 
which  has  received  a  great  deal  of  attention  in  recent  years  ( [ A84-1] ,  [A84-2],  [A84-3], 
[BF81],  [BL84],  [F80],  [F83],  [L77],  [L84J,  [Z70],  and  elsewhere).  The  case  of  C1  cubics 
(piecewise  cubic  polynomials  with  one  continuous  derivative)  has  received  a  good  deal  of 
this  attention.  In  this  paper,  a  method  of  cheaply  constructing  such  interpolants  is  given. 
The  idea  is  to  take  any  old  piecewise  cubic  interpolant  (which  need  not  be  C1)  and  compute 
a  perturbation  of  it  which  satisfies  the  desired  smoothness  conditions. 

The  method  outlined  in  this  paper  is  based  on  the  B-form  representation  of  bivariate 
polynomials,  as  given  in  [B86]  and  [F80].  Readers  not  familiar  with  the  B-form  should  read 
[B86]  for  a  detailed  description  of  it  in  its  most  general  form.  In  any  case,  the  features  of 
the  bivariate  B-form  required  to  establish  the  schemes  in  this  paper  will  be  reviewed.  See 
[B86j  for  proofs  of  these  facts. 

Consider  the  points  u,  v,  w  6  R2.  An  arbitrary  point  x  €  R2  can  be  expressed  as 

an  affine  combination  of  u,  v,  and  w,  i.e.  x  =  au(x)u  a.,(x)r  +  au,(x)u;,  where  au(x)  + 

o,,(x)  -f  Ou,(i)  -  1.  The  numbers  au(x).  a„(x).  and  au.(x)  are  called  the  barycentric 

coordinates  of  x  with  respect  to  u,  v,  and  w ,  and  they  are,  in  fact,  affine  functions  of  their 
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argument.  Denote  by  a  the  vector  of  barycentric  coordinates,  and  consider  the  mult  i  udex 
i  £  Z\  with  iu  +  iv  +  iw  =  n.  It  is  clear  that 


a 


U  V  w 

*  f#  ft  t 

1  U**V*U»’ 


(1) 


is  a  bivariate  polynomial  of  degree  at  most  n  (since  each  component  of  a  is  affine).  The 
(n^2)  valid  choices  for  »  give  rise  to  the  same  number  of  linearly  independent  bivariate 
polynomials  of  degree  <  n.  Since  the  dimension  or  the  space  of  bivariate  polynomials  of 
degree  <  n  is  also  (n*2),  the  functions  a*  must  form  a  basis  for  that  space.  Thus,  any 
polynomial  p  of  degree  <  n  can  be  written  in  the  form 


p  - 

i 

for  certain  coefficients  e,.  This  is  the  B-form  of  a  bivariate  polynomial. 

This  form  is  so  useful  because  it  has  a  nice  geometric  interpretation.  Consider  the 
triangle  given  as  the  convex  hull  of  u,  v,  and  w.  Then  all  the  valid  choices  of  t  can  be 
viewed  as  occupying  locations  on  a  certain  mesh  over  the  triangle.  This  mesh  is  the  one 
generated  by  considering  all  of  the  points  in  the  triangle  with  barycentric  coordinates  i  n, 
for  each  of  the  valid  choices  of  *.  For  the  case  n  =  3,  the  locations  are  given  in  Figure 
1.  For  example,  the  point  corresponding  to  u  itself  has  iu  =  3,  iv  =  tw  =  0.  The  point 
nearest  u  along  the  edge  from  u  to  v  has  *u  =  2,  iv  =  1,  and  iw  =  0.  The  point  in  the 
center  has  iu  =  iv  —  xw  =  1. 

This  form  ha^  the  nice  property  that  the  values  of  the  coefficients  say  something  about 
the  local  behavior  of  p.  For  example,  if  iu  =  n,  with  iv  =  iu.  =  0,  then  c,  =  p(u),  i-e. 
the  coefficient  associated  the  each  of  the  vertices  is  just  the  value  of  the  polynomial  at 


that  vertex.  Moreover,  if  the  restriction  of  p  to  the  line  determined  by  one  of  the  edges 
is  considered,  then  that  univariate  polynomial  is  uniquely  determined  by  the  coefficients 
associated  with  that  edge.  These  properties  make  it  possible  to  discuss  representations  of 
piecewise  polynomial  functions  using  a  generalization  of  the  B-form,  namely  the  B-net. 


Figure  1 

Consider  more  generally  4  points  in  R2,  say  u,  v,  w,  and  z,  whose  convex  hull  is  a 
quadrilateral.  Suppose  that  this  quadrilateral  is  divided  into  two  triangles,  one  of  which 
is  the  convex  hull  of  u,  v,  and  w,  and  the  other  is  the  convex  hull  of  v,  tv,  and  z.  Let 
the  quadrilateral  be  the  domain  for  some  piecewise  polynomial  function  of  degree  <  n 
with  two  pieces,  each  of  which  has  as  its  domain  one  of  the  triangles.  It  is  clear  that  each 
polynomial  may  be  described  in  B-form  on  each  triangle.  If,  as  is  generally  the  case,  the 
pp  function  is  continuous,  then  the  univariate  polynomials  which  are  the  restrictions  of 
each  polynomial  to  the  common  edge,  in  this  case  given  by  v  and  w ,  must  be  the  same. 
Since  these  polynomials  are  determined  by  the  B-form  coefficients  along  those  edges,  and 
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the  polynomials  must  be  the  same,  the  coefficients  must  also  be  the  same.  The  B-form 
representation  of  continuous  pp  functions  is  called  the  B-net.  The  B-net  for  the  above 
example,  with  n  =  3,  is  given  in  Figure  2. 


w 


Figure  2 

In  this  paper,  the  goal  is  to  construct  C1  interpolants  with  cubic  pp  functions.  Thus, 
simple  continuity  of  pp  functions  is  not  enough.  Fortunately,  the  conditions  for  C1  smooth¬ 
ness  in  terms  of  the  B-net  are  both  simple  and  elegant.  In  Figure  2  there  are  three  quadri¬ 
laterals  marked  with  a  dashed  line,  each  of  which  is  similar  to  the  domain.  For  the  pp 
function  given  to  be  C1,  it  is  necessary  and  sufficient  that  for  each  small  quadrilateral,  the 
four  B-net  coefficients  associated  with  its  vertices  satisfy  a  certain  linear  relationship:  If 
the  four  coefficients  are  viewed  as  function  values  at  the  points  u,  v ,  w,  and  z,  related  by 
the  similarity  of  the  small  quadrilaterals  to  the  big  one,  then  the  graph  of  the  four  points 


must  be  planar.  In  general,  n  of  these  conditions  will  be  required  for  pp  functions  of  degree 


<  n. 

To  turn  now  to  the  main  topic  of  the  paper,  consider  any  collection  of  points  in  R2 
at  which  the  values  of  some  function  which  is  to  be  interpolated  are  given.  Assume  that 
a  triangulation  of  these  points  is  also  given.  If  not,  one  can  always  be  computed  by  the 
method  outlined  in  [GS78J.  The  so-called  Delauney  triangulation  produced  by  this  method 
is  viewed  by  many  as  being  a  good  one  for  the  purpose  of  interpolation,  and  it  has  been 
shown  ( [Si78] )  to  be  the  same  triangulation  as  that  produced  by  a  very  different  method 
in  [L72j.  In  any  case,  if  there  are  d  points  in  the  interior  of  this  triangulation  and  e  on  the 
boundary,  then  there  will  be  2d  e  -  2  triangles  in  the  triangulation  [BL84j.  From  this 
fact,  it  follows  that  there  are  3d  +  2e  —  3  edges  in  the  triangulation,  of  which  3d  +  e  -  3 
are  interior  edges. 

Given  this  collection  of  points,  consider  piecewise  cubic  polynomial  functions  defined 
over  the  triangulation,  i.e.  functions  which  are  locally  cubic  polynomials  in  each  triangle. 
The  B-net  for  these  functions  consists  of  9d  +  6e  -  8  coefficients.  However,  not  all  of 
these  coefficients  are  variables  in  the  situation  considered  here.  Only  those  functions 
which  actually  interpolate  need  to  be  considered.  This  uniquely  determines  the  coefficients 
located  at  the  vertices,  leaving  only  8d  +  5e  —  8  variables.  Making  the  interpolant  C1 
amounts  to  imposing  three  linear  conditions  at  each  interior  edge,  for  a  total  of  9d  +  3e  -  9 
equations.  Thus,  computing  a  bivariate  C1  cubic  pp  interpolant  for  scattered  data  boils 
down  to  solving  the  linear  system 

Ax  =  b  (2) 
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for  some  9d +  3e-9  by  8d-f  5e  —  8  matrix  A  and  an  8d+5e  —  8  vector  b ,  both  of  which  depend 
upon  the  location  of  the  vertices  and  the  structure  of  the  triangulation.  The  vector  6  also 
depends  upon  the  function  values  which  are  to  be  interpolated.  The  smoothness  conditions 
are  homogeneous  except  when  one  of  the  four  coefficients  is  known.  This  happens  exactly 
when  one  of  the  four  points  is  a  vertex  of  the  triangulation. 

A  great  deal  is  known  about  A  and  b.  It  is  known,  for  example,  that  A  is  extremely 
sparse,  having  either  three  or  four  non-zero  entries  in  each  row  and  no  more  than  three 
non-zero  entries  in  each  column.  Exactly  one  third  of  the  entries  in  b  are  zero,  and  these 
correspond  to  the  rows  in  A  which  have  exactly  four  non-zero  entries.  It  is  also  known 
([MS77],  [S79],  and  [S84])  that  the  row  rank  of  A  is  no  larger  than  Id  +  e  -  9,  which 
guarantees  that  the  linear  system  (2)  is  underdetermined.  It  also  means  that  (2)  might 
not  be  solvable.  As  of  this  writing,  the  solvability  of  (2)  remains  a  conjecture.  No  one 
has  been  able  to  prove  the  existence  of  a  solution  ot  to  produce  a  set  of  points  (and  a 
triangulation)  for  which  there  is  none.  For  the  remainder  of  this  paper,  the  existence  of  a 
solution  to  the  linear  system  (2)  will  be  assumed. 

In  order  to  construct  C 1  cubic  pp  interpolants  to  the  data,  some  way  of  eliminating  the 
extra  d+4e  + 1  degrees  of  freedom  which  arise  because  of  the  rank  deficiency  of  A  is  needed. 
One  way  of  doing  this  is  to  begin  by  ignoring  the  linear  system  (2)  and  to  construct  any  old 
interpolant.  For  example,  the  piecewise  linear  interpolant  may  be  considered.  It  is  very 
easy  to  construct,  since  the  value  of  each  of  the  cubic  B-net  coefficients  is  simply  the  value 
of  the  piecewise  linear  interpolant  at  that  point  ((B86J  and  [F80]).  In  any  case,  suppose 
that  some  interpolant  is  given,  and  let  x  be  the  vector  of  B-net  coefficients  corresponding 
to  it.  It  is  clear  that,  in  general,  Ax  ±  6,  since  A  and  6  were  never  considered  in  forming 


x.  However,  this  does  provide  a  way  of  eliminating  the  extra  degrees  of  freedom  inherent 
in  (2).  AC1  piecewise  cubic  interpolant  to  the  given  data  can  be  computed  by  solving 


(3) 


min  ||e|| 2 
subject  to  Ae  =  r, 

where  r  :=  b  -  Ax,  and  A  and  b  are  as  before.  If  e  solves  (3),  then  A(x  +  e)  =  b,  so  solving 
(3)  amounts  to  computing  a  perturbation  of  the  original  interpolant  in  order  to  make  it 

C*. 


The  problem  two  may  be  solved  by  many  different  techniques.  One  way,  for  example, 
would  be  to  compute  the  singular  value  decomposition  of  A.  While  this  method  has 
many  highly  desirable  numerical  properties,  it  is  inappropriate  here  because  it  fails  to 
take  advantage  of  the  sparsity  of  the  matrix  A.  The  same  is  true  of  Lemke’s  method  and 
other  pivotal  methods,  which,  even  if  they  make  use  of  sparse  matrix  techniques,  will  still 
depend  in  some  critical  way  on  the  ordering  of  the  rows  and  columns  of  A.  This  in  turn  will 
depend  upon  how  the  vertices  of  the  triangulation  (and  the  triangles  themselves)  happen 
to  be  ordered.  In  this  paper,  this  difficult  question  will  be  overlooked  by  using  an  iterative 
method,  namely  successive  overrelaxation  (SOR),  to  solve  the  problem. 

The  space  in  which  e  lies  can  be  decomposed  into  two  orthogonal  subspaces,  namely 
ker(A),  the  nullspace  of  A,  and  ran(Ar).  the  range  of  AT .  Any  solution  to  Ae  =  r  has  a 
unique  decomposition  e  =  e*  +  er,  where  e*  f  ker(A)  and  er  •£  ran(Ar).  In  order  to  obtain 
a  least  norm  solution  to  At  —  r,  i.e..  to  have  e  solve  (3).  ’.t  is  necessary  and  sufficient  to 
have  tk  =  0.  Thus,  if  e  :=  ATy,  then  solutions  to 


AATy  -  r 
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H) 


will  also  provide  solutions  to  (3). 


$ 


Let  L  be  the  matrix  the  same  size  as  AAT  whose  entries  are  the  same  as  AAT  below 
the  main  diagonal  and  zero  everywhere  else.  More  straightforwardly,  L  is  the  strict  lower 
triangular  part  of  AAT .  Let  D  be  the  diagonal  part  of  AAT,  so  that  L  +  D  +  LT  =  AAT . 
It  is  clear  that  D  has  strictly  positive  entries  on  its  main  diagonal,  since  each  of  these 


ins 

U 


entries  is  just  the  inner  product  of  a  row  of  A  with  itself.  Then  the  iterative  scheme 


yn  +  1  =yn  +  wD~l(r  -  Lyn+1  -  Dyn  -  LTyn) 


converges  for  any  choice  of  0  <  u  <  2.  Here  y”  is  the  vector  which  is  the  n-th  iterate  of  the 
scheme,  and  the  convergence  is  to  any  solution  of  the  system  (4).  Since  A  is  rank  deficient, 
it  is  clear  that  AAT ,  although  symmetric  and  positive  semi-definite,  may  nevertheless  be 
singular.  The  convergence  is  a  consequence  of  the  following  theorem,  due  to  Keller: 

Theorem  1  [K65j:  Let  5  be  a  symmetric  matrix  of  order  m  and  let  TV  be  a  non- 
singular  matrix  of  order  m  for  which 


P:=  Ni-  NT  -  S 


is  positive  definite.  Then  the  iterative  scheme 


(.V  S)x”  *  b 


converges  if  and  only  if  S  is  positive  semi-definite  and  b  ran (5).  The  convergence  is  to 
a  solution  of  Sx  =  b. 


i 


The  iterative  scheme  (5)  is  of  the  form  (6)  if  N  is  chosen  to  be  L  -f  ~D.  For  any  finite 

choice  of  u,  this  is  non-singular.  Then 

P  =  N  +  Nt  -  AAt 

=  (L  +  - D )  +  {L+-D)T  -{L  +  D+  Lt) 
uj  u> 

2  -  u)  _ 

=  - D, 

u) 

which  is  positive  definite  if  0  <  u  <  2.  Hence,  SOR  converges  to  some  solution  of  the 
linear  system  (4). 

As  it  turns  out,  it  doesn’t  matter  which  solution  to  (4)  is  obtained.  Any  solution  y 
has  two  orthogonal  components,  y *  £  ker(AAr)  and  yT  £  ran((AA  r)T)  =  ran(AAT),  so 
that  y  =  yk  +  yr-  However,  no  matter  what  yk  is,  ATyk  =  0,  since  ATyk  £  ran(Ar)  and 
ATyk  £  ker(A).  Thus,  c  =  ATy  depends  only  on  yr,  and  not  on  y*..  This  uniquely  solves 
the  system  (3). 

The  only  remaining  step  in  the  interpolation  process  is  to  take  e  and  add  it  to  x.  This 
will  provide  B-net  coefficients  which  satisfy  the  conditions  for  C1  smoothness,  interpolate 
the  given  data,  and  are,  in  the  2-norm  sense,  as  close  as  possible  to  those  in  the  original 
interpolant.  Since  the  polynomials  (l)  are  a  well-conditioned  basis  for  the  space  of  pp 
functions  considered  [H82],  it  is  reasonable  to  expect  that  the  resulting  interpolan,  will 
also,  in  some  sense,  be  as  close  as  possible  to  the  original  interpolant.  This  means  that 
if  the  original  interpolant  was  constructed  to  satisfy  certain  shape  requirements,  it  is  at 
least  conceivable  that  the  resulting  C1  interpolant  may  satisfy  those  same  r  juirements. 

Before  proceeding  to  the  numerical  experiments,  there  are  a  few  features  of  this  scheme 
that  are  worth  noting.  The  first  is  that  the  matrix  AAT  is  extremely  ,  arse.  Since  each 
row  of  A  has  at  most  four  non-zero  elements,  each  row  (and  column)  of  AAT  has  at  most 


eight  non-zero  elements  in  addition  to  the  entry  on  the  main  diagonal.  Because  it  is  also 
symmetric  only  half  ot  it  needs  to  be  stored  in  memory  in  any  computer  implementation 
of  the  scheme. 

It  is  also  worth  noting  that  the  assumed  solvability  of  the  linear  system  (2)  is  an 
essential  ingredient  in  Theorem  1.  Without  this,  the  iterative  scheme  cannot  converge. 
Since  this  is  the  only  condition  in  the  theorem  which  is  not  ironclad,  the  appearance  of  a 
problem  for  which  the  method  fails  to  converge  provides  a  candidate  for  a  counter-example 
to  the  conjecture.  In  all  experiments  performed  so  far,  non-convergence  has  never  been 
observed. 

The  numerical  experiments  which  have  been  performed  have  all  used  a  local  imple¬ 
mentation  of  the  scheme  (for  piecewise  linear  initial-guess  interpolants)  on  a  VAX-1 1/780 
running  VMS.  This  implementation  is  in  the  form  of  a  FORTRAN  subroutine  which  takes 
as  input  an  array  of  points  in  R2  and  an  array  of  corresponding  function  values.  It  begins 
by  calling  a  Delauney  triangulation  routine  to  determine  a  suitable  triangulation  of  the 
region  of  interpolation.  Next,  it  determines  the  piecewise  linear  interpolant  for  the  given 
data,  i.e.  the  vector  x.  After  forming  the  matrix  A  and  the  vector  b ,  the  matrix  AAT  and 
the  vector  r  =  6  —  Ax  are  computed.  SOR  is  then  used  to  solve  the  system  AAT y  =  r  to 
a  tolerance  specified  by  the  user  with  the  parameter  u;  also  specified  by  the  user.  Finally, 
the  vector  x  +  ATy  is  computed,  and  its  components  are  put  into  an  array  which  contains 
the  B-net  for  the  interpolant. 

One  of  the  most  important  features  for  any  interpolant  to  have  is  that  it  be  local, 
which  means  that  a  change  in  one  of  the  function  values  will  only  change  the  resulting 
interpolant  near  the  location  of  the  function  value.  This  property  has  been  investigated 


I 


for  a  number  of  different  data  sets  by  choosing  all  of  the  function  values  to  be  zero  except 


at  one  point  where  it  is  chosen  to  be  one.  The  graphs  of  the  resulting  interpolants  are  all 


very  similar  to  the  one  shown  in  Figure  3.  It  is  clear  from  this  that  the  interpolant  is  not 


local,  but  might  be  essentially  local,  meaning  that  the  effects  of  a  change  in  a  function 


value  decay  with  distance. 


If,  in  fact,  the  interpolant  is  essentially  local,  the  error  in  the  interpolant  should  be  at 


least  as  good  as  0(p2),  where  p  is  the  maximum  diameter  of  any  of  the  triangles  in  the 


triangulation.  Here,  diameter  means  the  length  of  the  longest  segment  which  is  contained 


within  a  triangle.  The  error  estimate  should  be  this  good  because  the  interpolation  scheme 


used  reproduces  linear  polynomials.  In  other  words,  if  the  given  data  describes  a  linear 


function,  then  the  interpolant  will  be  precisely  that  function. 


For  the  bivariate  function  f(x i,x2)  =  x\  +  x\  -  2xjx2  +  x^  +  2z2  +  3  defined  over 


the  square  [0,  l]2,  interpolants  on  successively  finer  triangulations  have  been  computed.  In 


each  case,  the  interpolation  points  are  located  at  the  vertices  of  an  m  x  m  square  mesh 


placed  over  [0,  l]2.  The  following  table  summarizes  the  results.  The  column  labelled  k  is 


the  observed  rate  of  convergence  between  each  entry  and  the  subsequent  one.  Here  k  is  the 


exponent  if  the  error  is  assumed  to  behave  like  0(pk).  The  errors  listed  are  the  maximum 


absolute  differences  between  the  function  and  its  interpolant  which  are  encountered. 


m  Error 


3  0.14779282  1.939 


4  0.06732988  1.989 


5  0.03799295  1.994 


mmmssssmam 


•tv  I* 


6 


0.02434874 


1.994 


7 

0.01692724 

1.998 

8 

0.01243925 

1.998 

9 

0.00952625 

2.015 

10 

0.00751376 

2.135 

12 

0.00489593 

1.831 

14 

0.00360560 

1.993 

16 

0.00271082 

The  observed  rate  of  convergence  is  very  nearly  the  0(p2)  which  had  been  hypothesized. 
It  would  have  been  shocking,  in  fact,  to  have  obtained  anything  better,  since  the  function 
is  quadratic  and  is  clearly  not  being  reproduced  by  the  scheme! 

Also  deserving  of  mention  is  the  fact  that  the  m  =  16  case,  which  is  interpolation 
of  256  function  values  with  450  (perhaps)  different  bivariate  cubic  polynomials,  required 
less  than  321.94  seconds  to  compute.  Given  that  this  involved  the  triangulation  of  the  256 
points  (a  needless  waste,  since  the  mesh  is  regular  anyway),  the  construction  and  solution 
of  a  1935  x  1935  linear  system  of  equations,  and  the  evaluation  of  the  interpoiant  on  a 
51  x  51  grid  (2601  interpoiant  evaluations),  this  seems  quite  reasonable.  Certainly  much 
larger  problems  are  possible  even  on  the  VAX. 

I  Of  course,  the  method  works  well  independent  of  the  regularity  of  the  triangulation. 

|  Consider  the  function  g{xx,  i2)  =  (ij  -  xf)(x2  -  x^e31*  _7l=.  A  contour  map  of  this 

|  function  over  the  square  [0, 1  ] 2  is  given  in  Figure  4.  For  this  function,  its  interpoiant  has 

|  been  computed  over  the  triangulation  given  in  Figure  5.  This  triangulation  is  the  Delauney 

I 

1 


Figure  5 


Figure  7 
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